COMMENT ⊗   VALID 00003 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	(*PROGRAM HEADER PAGE*)
C00005 00003	PROGRAM FFT(OUTPUT)
C00013 ENDMK
C⊗;
(*PROGRAM HEADER PAGE*)

(*PAS10 OPTIONS*) (* D+,R32*)

(*	            					     DEFAULT 

D+	DEBUG AND POSTMORTEM DUMP				-
E+	EXTERNAL CALLS TO LEVEL 1 PROCEDURES ALLOWED		-
Fn	FILE OPTION						1
I+	FORTRAN I/O IN EXTERNAL FORTRAN SUBROUTINES		-
L+	OBJECT LISTING						-
Rn 	SIZE OF LOW-SEGMENT 			        (SEE PAS10 MANUAL)
Sn   	MAX INSTRUCTIONS PER STATEMENT			       1000
T+	RUNTIME CHECK						+
U+	72 COLUMN FORMAT					-
Xn	HIGHEST REGISTER FOR PARAMETERS				6
*)

(*SLAC PCPASC OPTIONS*) (* B+,D+,M-*)

(*                   					     DEFAULT

A+	GENERATE 370 OBJECT MODULE				-
A-	GENERATE 370 ASSEMBLY MODULE
B+	BOUNDS CHECKING, BUT ALLOW 'BIG' CHARACTERS		-
C+	EMIT PCODE						+
D+	RUNTIME CHECKING OF POINTER, INDEX, SUBRANGE VALUES	-
E+	FILE IS IN EBCDIC CHARACTER SET				-
F+	SAVE FPR'S ON PROCEDURE/FUNCTION ENTRY			+
K+	ENABLE STATEMENT EXECUTION COUNTING			-
L+	LIST SOURCE PROGRAM					+
M+	72 COLUMN FORMAT					+
P+	DOUBLE-WORD BOUNDARY ALIGNMENT				-
S+	SAVE GPR'S ON PROCEDURE/FUNCTION ENTRY			+
T+	PRINT SYMBOL TABLES (FOR POST-PROCESSOR)		-
U+	GET STATISTICS?? 2ND PARAMETER TO PCODE BGN INSTR.	-
V+	?? 3RD PCODE BGN INSTRUCTION PARAMETER			-
X+	USE ACTUAL PROCEDURE NAMES FOR EXTERNAL REFERENCES	-
X-	GENERATE UNIQUE 8-CHAR NAMES FOR EXTERNAL REFERENCES
*)

(*S1 PCPASC OPTION DIFFERENCES*) (*$A+,B+,D+,M-*)

(*							     DEFAULT
	
A+	GENERATE S1 ASSEMBLY MODULE				-
A-	GENERATE S1 OBJECT MODULE
*)

(* SLAC/PDP-10 TRANSPORT DEPENDENCIES FLAGGED WITH "XSL10" *)
(* PDP-10/S-1 TRANSPORT DEPENDENCIES FLAGGED WITH "X10S1" *)
PROGRAM FFT(OUTPUT);

CONST
M=	6;
N=	64;
NLINES= 64;
NPRD=	4;
NCOL=	35;
PI=	3.1415926536;

TYPE
COMPLEX=	RECORD RL: REAL; IM: REAL END;

VAR
XARG:		REAL;
SIG:		ARRAY [1 .. N] OF COMPLEX;
TRANS:		ARRAY [1 .. N] OF COMPLEX;

(**********************************************************************)

PROCEDURE CADD(VAR RESULT: COMPLEX;
	       X, Y: COMPLEX);
BEGIN
RESULT.RL:=X.RL+Y.RL;
RESULT.IM:=X.IM+Y.IM;
END;

(**********************************************************************)

PROCEDURE CSUB(VAR RESULT: COMPLEX;
	       X, Y: COMPLEX);
BEGIN
RESULT.RL:=X.RL-Y.RL;
RESULT.IM:=X.IM-Y.IM;
END;

(**********************************************************************)

PROCEDURE CMULT(VAR RESULT: COMPLEX;
	        X, Y: COMPLEX);
BEGIN
RESULT.RL:=X.RL*Y.RL-X.IM*Y.IM;
RESULT.IM:=X.RL*Y.IM+X.IM*Y.RL;
END;

(**********************************************************************)

FUNCTION PWR2(EXP: INTEGER): INTEGER;

VAR
I,P :  INTEGER;

BEGIN
P:=1;
FOR I:=1 TO EXP DO P:=P*2;
PWR2:=P;
END;
 
(**********************************************************************)

FUNCTION LLWBOTH: REAL;
VAR
X,SIGNREV,XSQ :  REAL;

BEGIN

X := XARG;
SIGNREV := 1.0;
IF X < 0.0 THEN BEGIN
  X := -X;
  SIGNREV := -SIGNREV
END;

X := X - 2*PI*TRUNC(X/(2*PI));
CASE TRUNC(X/(PI/2)) OF
 0: ;
 1: X := PI - X;
 2: BEGIN X := X - PI;  SIGNREV := -SIGNREV END;
 3: BEGIN X := 2*PI - X;  SIGNREV := -SIGNREV END;
 4: X := X - 2*PI
END (*case*);

XSQ := X*X;
LLWBOTH := SIGNREV*X*(1-XSQ*(1/6-XSQ*(1/120-XSQ*(1/5040))));
(*  LLWBOTH := SIGNREV*X*(1-XSQ*(IFACT3-XSQ*(IFACT5-XSQ*(IFACT7)))); *)

END (*LLWBOTH*);

(**********************************************************************)

FUNCTION LLWCOS(X: REAL): REAL;

BEGIN
XARG:=X+PI/2;
LLWCOS:=LLWBOTH
END;

(**********************************************************************)

FUNCTION LLWSIN(X: REAL): REAL;

BEGIN
XARG:=X;
LLWSIN:=LLWBOTH
END;

(**********************************************************************)

FUNCTION LLWSQRT(X: REAL): REAL;

CONST
THRESH=	0.00000001;

VAR
APPROX :  REAL;
LAST←APPROX :  REAL;

BEGIN

APPROX:=X;
REPEAT
  LAST←APPROX:=APPROX;
  APPROX:=(APPROX+X/APPROX)/2;
UNTIL ABS(APPROX-LAST←APPROX)<THRESH;
LLWSQRT:=APPROX;

END;

(**********************************************************************)

PROCEDURE INITSIG;

VAR
I :  INTEGER;

BEGIN
  FOR I:=1 TO N DO BEGIN
    SIG[I].RL:=LLWSIN((I-1)*2*PI*NPRD/N);
    SIG[I].IM:=0;
    TRANS[I]:=SIG[I];
  END;
END;

(**********************************************************************)

PROCEDURE WRT←GRAPHS;

VAR
J,X←INDX,I,NSPACE :  INTEGER;
MAX←SIG,
MIN←SIG,
SIG←SCALE,
MAX←TRANS,
MIN←TRANS,
TRANS←SCALE :  REAL;

BEGIN

FOR I:=1 TO N DO
  TRANS[I].RL:=LLWSQRT(TRANS[I].RL*TRANS[I].RL+TRANS[I].IM*TRANS[I].IM);

MAX←SIG:=SIG[1].RL;
MIN←SIG:=SIG[1].RL;
MAX←TRANS:=TRANS[1].RL;
MIN←TRANS:=TRANS[1].RL;

FOR I:=2 TO N DO BEGIN
  IF SIG[I].RL>MAX←SIG THEN MAX←SIG:=SIG[I].RL;
  IF SIG[I].RL<MIN←SIG THEN MIN←SIG:=SIG[I].RL;
  IF TRANS[I].RL>MAX←TRANS THEN MAX←TRANS:=TRANS[I].RL;
  IF TRANS[I].RL<MIN←TRANS THEN MIN←TRANS:=TRANS[I].RL;
END;

IF (MAX←SIG-MIN←SIG)<1.0E-20
THEN SIG←SCALE:=1.0
ELSE SIG←SCALE:=NCOL/(MAX←SIG-MIN←SIG);

IF (MAX←TRANS-MIN←TRANS)<1.0E-20
THEN TRANS←SCALE:=1.0
ELSE TRANS←SCALE:=NCOL/(MAX←TRANS-MIN←TRANS);

FOR I:=1 TO NLINES DO BEGIN
  X←INDX:=TRUNC(N*I/NLINES);
  
  NSPACE:=TRUNC(SIG←SCALE*(SIG[X←INDX].RL-MIN←SIG));
  FOR J:=1 TO NSPACE DO WRITE(' ');
  WRITE('*');
  FOR J:=NSPACE-2 TO NCOL DO WRITE(' ');

  NSPACE:=TRUNC(TRANS←SCALE*(TRANS[X←INDX].RL-MIN←TRANS))+3;
  FOR J:=1 TO NSPACE DO WRITE(' ');
  WRITE('*');

  WRITELN();
END;
  
END;

(**********************************************************************)

PROCEDURE XFORM;

VAR
I,J,K,L,IP,LE,LE1 : INTEGER;
U,W,T :  COMPLEX;

BEGIN

J:=1;
FOR I:=1 TO N-1 DO BEGIN
  IF I<J THEN BEGIN
    T:=TRANS[J];
    TRANS[J]:=TRANS[I];
    TRANS[I]:=T;
  END;
  K:=N DIV 2;
  WHILE K<J DO BEGIN
    J:=J-K;
    K:=K DIV 2;
  END;
  J:=J+K;
END;

FOR L:=1 TO M DO BEGIN
  LE:=PWR2(L);
  LE1:=LE DIV 2;

  U.RL:=1.0;  
  U.IM:=0.0;

  W.RL:=LLWCOS(PI/LE1);
  W.IM:=LLWSIN(PI/LE1);

  FOR J:=1 TO LE1 DO BEGIN
    I:=J;
    REPEAT
      IP:=I+LE1;
      CMULT(T,TRANS[IP],U);
      CSUB(TRANS[IP],TRANS[I],T);
      CADD(TRANS[I],TRANS[I],T);
      I:=I+LE;
    UNTIL I>N;
    CMULT(U,U,W);
  END;
END;

END;

(**********************************************************************)

BEGIN

INITSIG;
XFORM;
WRT←GRAPHS;

END.

(**********************************************************************)